Coordinate systems and maps

Goals:

  1. In the first exercise, you use the Grammar of Graphics to fix a map with tmap() by modifying the provided code. The libraries, data, and initial code is provided below.

  2. In the second exercise, you inspect and assess the compatibility of the provided spatial datasets

  3. In the third exercise, you

  • extrapolate a missing CRS, and
  • define the missing projection in R
  1. In the fourth (and optional) exercise, practice how to
  • reconcile CRS of your spatial data
  • reproject raster and vectors to a shared SRS
  1. In the fifth exercise, improve your a map of Aarhus from Week 2 following the best map design principles.
  • plot layers together
  • limit your data into an area of interest
  1. Create a mirror pair of honest and lying maps of Danish 2025 kommunal election results

Required R libraries

We will use the sf, raster/terra, and tmap packages and then a few oldies.

library(sf)
library(raster)
library(tmap)
library(googlesheets4)
library(tidyverse)

Exercise 1 - Fix a map of Denmark

In this exercise you will learn to make a map with tmap library, by adding spatial data layers and modifying the arguments that specify their rendering

Data sets

We will use two data sets: elevation and counties. Both originate from the UC Davis global databaset and you can redownload them with getData(). (see Week02, Exercise 12) The first one is an elevation raster object for Denmark, and the second one is a spatial dataframe object with polygons representing the 99 Danish municipalities.

Existing code

Here is the code to create a new map of Denmark. Your role is to improve this map based on the suggestions below.

# get elevation
elevation <- raster("../data/DNK_msk_alt.grd")

# get counties
counties <- readRDS("../data/gadm36_DNK_2_sp.rds")


# starting a static map in tmap v3
tm_shape(elevation) + tm_raster(title = "elev", style = "cont", palette = "BuGn") +
    tm_shape(counties) + tm_borders(col = "red", lwd = 3) + tm_scale_bar(breaks = c(0,
    100, 200), text.size = 1) + tm_compass(position = c("LEFT", "center"), type = "rose",
    size = 2) + tm_credits(text = "A. Sobotkova, 2024") + tm_layout(main.title = "My terrible map",
    bg.color = "orange", inner.margins = c(0, 0, 0, 0))

# starting a static map in tmap v4
tm_shape(elevation) + tm_raster(col.legend = tm_legend(title = "elev"), col.scale = tm_scale_continuous(values = "brewer.bu_gn",
    midpoint = 0)) + tm_shape(counties) + tm_borders(col = "red", lwd = 3) + tm_scalebar(breaks = c(0,
    100, 200), text.size = 1) + tm_compass(position = c("LEFT", "center"), type = "rose",
    size = 2) + tm_credits(text = "A. Sobotkova, 2025") + tm_title("My terrible map") +
    tm_layout(bg.color = "orange", inner.margins = c(0, 0, 0, 0))

# install.packages(c('shinyjs', 'colorblindcheck')) cols4all::c4a_gui()

Tasks

  1. Change the map title from “My terrible map” to “Denmark’s municipalities”.
  2. Update the map credits with your own name and today’s date.
  3. Change the color palette to yellow to brown sequence (“-RdYlGn” in v3 and equivalent in v4). (You can also try other palettes from http://colorbrewer2.org/ or cols4all::c4a_gui())
  4. Put the north arrow in the top right corner of the map or some better position.
  5. Improve the legend title by adding the used units (masl).
  6. Increase the number of breaks in the scale bar and move it to top right.
  7. Change the borders’ color of the municipalities to black. Decrease the line width.
  8. Change the background color to some reasonable color of your choice.

Your solution

# /Start Code/

# /End Code/

Exercise 2 - Create and inspect spatial data

We will use three data sets: elevation, places and nitrates, explore them, and make a map. The first one is an elevation raster object for Denmark you already used above. The second dataset is a list of places you created in week 1. The third is a geochemical dataset from landbrug.dk.

You have already read in the raster, and so now the focus is on the vector data. You need to create these from googlesheet/csv before you can inspect them. Follow the instructions and answer the questions related.

Preparation: Create vector data from a csv with LatLong and WKT

Load the two datasets:

  • read in the places googlesheet using the read_sheet() function (grab it from W01 exercise)
  • use read_csv() to grab nitrates.csv file and use slice() function to grab the first 5000 records upon reading it in to reduce its size.
# Load your DK places googlesheet from W01
gs4_deauth() # if the Gdrive authentication is not working for you

places <- read_sheet("https://docs.google.com/spreadsheets/d/1PlxsPElZML8LZKyXbqdAYeQCDIvDps2McZx1cTVWSzI/edit#gid=1817942479",_______)


# Load the first 5000 rows in .csv-file and save it as "nitrates":
nitrates <- read_csv("../data/nitrate.csv") %>% 
  slice(1:5000)

Convert the tabular data into simple feature using their geometry columns

Tabular data contain spatial data, that you can use to convert these tibbles into simple features! Note the lat/long columns in places and wkt column in nitrates. Inspect these columns and convert the objects to simple feature using the st_as_sf() function using the coords and the wkt arguments respectively!

  • take places, filter missing coordinates away and then use the st_as_sf() function with the coords argument to cast into a simple feature places_sf
  • for nitrates, use the st_as_sf() function with the wkt argument to convert to simple feature nitrates_sf.
  • use plot(places_sf$geometry) and plot(nitrates_sf$WKT) to view each simple feature in space.
places_sf <- places %>% 
  filter(!is.na(_________)) %>% 
  st_as_sf(coords = ____(__________,________))

nitrates_sf <- _______(nitrates, wkt = __________)

Instructions and questions

Type answers to the questions as code comments next to or under the code used

  • Display the nitrates_sf and places_sf objects and view their structure. What can you say about the content of each file?

    • What type of data does it store?
    • How many attributes does it contain?
    • What is its geometry? (points, lines, polygons, surfaces, network)
    • What kind of coordinate system does it use?

  • Display the elevation object and view its structure. What can you say about the content of this file?

    • What type of data does it store?
    • What is the coordinate system used?
    • How many attributes does it contain?
    • How many dimensions does it have?
    • What is the data resolution?

  • Can you plot the files on top of one another? Why?

Your solution and answers

# /Start Code / #


# /End Code/ #

Exercise 3 - Define coordinate systems

A machine-readable CRS is prerequisite to spatial visualisation and analysis, but not all spatial datasets that float around have its CRS defined. Here you learn to define/specify the coordinate systems for datasets where the CRS is NA/missing/unrecognized by R. You can fix this situation if YOU know what the CRS should be.

Look at the coordinates in places and nitratesagain. While both refer to places in Denmark, their formats are different, signalling different CRS. Can you guess what they are?

  • places contain lat/long columns which you have created on the basis of GoogleMaps which uses Web Mercator, or WGS84 CRS. This translates into EPSG 4326.
  • nitrates contain WKT column with planar coordinates provided by landbrug.dk server. The source server specifies that all data is ETRS89 UTM32N, which corresponds to EPSG 25832.

Instructions:

  • define CRS for each spatial object with st_set_crs() using the CRS they are compatible with.
  • write the CRS into the name of these defined object in R and write the CRS into the name!

Your solution and answers

# /Start Code / #

places4326 <- places_sf %>% _______(_________)

nitrates25832 <- nitrates_sf %>% _________(________)

places4326
nitrates25832
# /End Code/ #

Wonderful, your vector data now have defined spatial metadata. In the next step, reconcile these two CRSs so you can display them in the same map.

Exercise 4 - Reconce coordinate systems

Now that you know that coordinate systems differ, make them compatible. You have two final CRS choices. Convert your objects to one or another CRS so that you can display the nitrate samples and DK places over the elevation map of Denmark.

Remember that you use st_transform() to reproject vector data, and projectRaster() on raster data to change their CRS.

Instructions

  • Option 1: Transform the nitrates dataset into the coordinate reference system used in the places object.

    • Create a new object nitrates with the places crs. You can label it nitrates_#### writing the EPSG out for easy differentiation.
    • Visualize the results using the plot() function.

  • Option 2: Reproject the elevation and places data into the coordinate reference system used in the nitrates object.

    • Create a new object elevation_#### with the nitrates crs.
    • Visualize the results (elevation_####`` together withnitrates`)

Your solution and answers

# /Start Code / #


# /End Code/ #

Exercise 5 - Improve your map of Aarhus from Week 2

Make your Week 2 map publication-ready! Download data about DK available online (if needed). Make a clean and legible map of Aarhus and a focus area (see below) with tmap library following the Exercise no.1 above and the guide here or here

Additional challenge:

  • Draw a bounding box (hint st_make_grid() to generate it) or a buffer around the district you live in or some other area of interest, and generate a more detailed map of this area, while leaving Aarhus + Denmark’s boundary in an inset for orientation. Make sure to add title and an explanation for your choice. Include a legend for any environmental background data.
  • Choose a reasonable classification if you are visualising environmental or quantitative data of any sort.
  • Figure out how you could print this result as .png ( explore ggsave) so the ratio and position of the inset is preserved

Your solution

# /Start Code/ #

# /End Code/ #

Exercise 6 - Challenge: Make a map that is true and that lies

You are a political consultant. Your client wants a map of Danish election results (2025) that makes their party look strong. Choose your party (e.g. DF) and create:

  1. An honest map: Follows cartographic best practices
  2. A misleading map: Technically correct but emphasizes your client’s strength (through color, classification, or geographic focus)
  3. A 300-word memo: Explains how the misleading map manipulates perception without lying

Get data

# get municipalities
muns <- readRDS("../data/gadm36_DNK_2_sp.rds") %>%
    st_as_sf()

# get votes
library(googlesheets4)
gs4_deauth()
votes <- read_sheet("https://docs.google.com/spreadsheets/d/1R77xvSEdEF5UMF6u1POlHVyYHbB1zFkNqEYTkwD7GPw/edit?")

# align municipality names
`%nin%` <- Negate(`%in%`)
votes$Division %in% muns$NAME_2
 [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
[25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
[49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[73]  TRUE  TRUE  TRUE
 [ reached getOption("max.print") -- omitted 23 entries ]
votes$Division[votes$Division %nin% muns$NAME_2]
[1] "Aarhus"        "Nordfyn"       "Høje-Taastrup" "Copenhagen"   
muns$NAME_2[31] <- "Aarhus"
muns$NAME_2[21] <- "Høje-Taastrup"
muns$NAME_2[25] <- "Copenhagen"

# join up

mun_vote <- muns %>%
    select(NAME_2) %>%
    left_join(votes, by = c(NAME_2 = "Division"))
mun_vote
Simple feature collection with 99 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 8.076389 ymin: 54.55903 xmax: 15.19306 ymax: 57.75153
Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:
       NAME_2    A    B    C    F   I   M    O    V   Æ    Ø   Å Others
1 Albertslund 25.3 11.3  7.4 16.1 3.8 1.2 11.9  3.3 1.5 14.4 2.1    1.8
2     Allerød 14.2  9.5 21.5 12.4 6.1 2.0  8.0 17.7 1.1  5.9 0.8    0.8
3    Ballerup 33.9  7.9  8.9 10.6 6.1 2.4 13.7  5.6 1.6  7.2 1.1    1.1
4    Bornholm 46.0  2.6  5.5  5.6 1.9 0.6 16.2  7.6 2.7  7.1 0.5    3.6
5     Brøndby 28.3 10.2  8.3 10.3 4.2 2.1 17.0  4.2 2.0 10.6 1.6    1.4
                        geometry
1 MULTIPOLYGON (((12.37129 55...
2 MULTIPOLYGON (((12.20744 55...
3 MULTIPOLYGON (((12.42482 55...
4 MULTIPOLYGON (((14.91719 55...
5 MULTIPOLYGON (((12.44434 55...
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

Look at the data

Is everything in place, any municipalities missing?

# check out the voting data
library(mapview)
mapview(mun_vote, zcol = "A")

Choose a map that is honest and one that lies

a <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "pretty", n = 5)
a

# 'Jenks' style further smooths over the gaps
b <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "jenks", n = 6)
b

# quantile style divides into 5 even groups
c <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "quantile", n = 5)
c

# Equal interval style divides the distribution into even groups
d <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "equal", n = 5)
d

# Write maps above to objects and plot them side by side with tmap_arrange()
# for better comparison

tmap_arrange(a, b, c, d)

Refs

Tennekes, Martijn. 2019. Tmap: Thematic Maps. https://CRAN.R-project.org/package=tmap.